setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM")
library(openxlsx)
sci=read.csv("E:/0 公共数据库差异情况/2018 Science supp/aat8127_Table_S1a_DGE.csv",header=T)

sci_scz=sci[sci$SCZ.fdr<0.1,]
sci_bd=sci[sci$BD.fdr<0.1,]
sci_scz=data.frame(symbol=sci_scz$gene_name,SCZ.2018sci.fdr=sci_scz$SCZ.fdr)
sci_bd=data.frame(symbol=sci_bd$gene_name,BD.2018sci.fdr=sci_bd$BD.fdr)
sci_scz$group.2018sci.scz="SCZ"
sci_bd$group.2018sci.BD="BD"

cmc.nosva=read.xlsx("E:/0 公共数据库差异情况/CMC/DEG_nosva.xlsx",sheet=1)
cmc.nosva=cmc.nosva[cmc.nosva$adj.P.Val<0.1,]
cmc.nosva=cmc.nosva[!cmc.nosva$MAPPED_genes==".",]
cmc.nosva=data.frame(symbol=cmc.nosva$MAPPED_genes,cmc.nosva.fdr=cmc.nosva$adj.P.Val)
cmc.nosva$cmc.nosva.group="CMC.nosva"

cmc.sva=read.xlsx("E:/0 公共数据库差异情况/CMC/DEG_sva.xlsx",sheet=1)
cmc.sva=cmc.sva[cmc.sva$adj.P.Val<0.1,]
cmc.sva=cmc.sva[!cmc.sva$MAPPED_genes==".",]
cmc.sva=data.frame(symbol=cmc.sva$MAPPED_genes,cmc.sva.fdr=cmc.sva$adj.P.Val)
cmc.sva$cmc.sva.group="CMC.sva"

file=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201222/807ASH.add.eQTL.GWAS.csv",header = T)
file=tidyr::separate_rows(file,Gene.refGene,sep=";")
file=merge(file,sci_scz,by.x="Gene.refGene",by.y = "symbol",all.x = T)
file=merge(file,sci_bd,by.x="Gene.refGene",by.y = "symbol",all.x = T)
file=merge(file,cmc.nosva,by.x="Gene.refGene",by.y = "symbol",all.x = T)
#file=merge(file,cmc.sva,by.x="Gene.refGene",by.y = "symbol",all.x = T)
file$DEG.validate.num=rowSums(file[,c(seq(189,194,2))]<0.1,na.rm = TRUE)
write.csv(file,"E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201216/at.least.two.ASH.in.DC.BF.10.add.eQTL.GWAS.add.DEG.csv",quote=F,row.names = F)

file1=file[!is.na(file$Other),]###eQTL.data
file1=file1[file1$DEG.validate.num>0,]
file1=file1[file1$DC.FDR.sig>4,]

length(unique(file1$unitID))
sel=which(file1$DC.pattern.valid.CC=="FALSE"|file1$DC.pattern.valid.HC=="FALSE")
file1=file1[-sel,]
